library(tidyverse)
library(Seurat)
library(cowplot)
library(ComplexHeatmap)
library(circlize)
library(GeneOverlap)
library(gprofiler2)
library(ggrepel)
library(ggplot2)
library(muscat)
library(purrr)
library(limma)
library(scran)
── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ── ✔ dplyr 1.1.0 ✔ readr 2.1.4 ✔ forcats 1.0.0 ✔ stringr 1.5.0 ✔ ggplot2 3.4.1 ✔ tibble 3.1.8 ✔ lubridate 1.9.2 ✔ tidyr 1.3.0 ✔ purrr 1.0.1 ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ── ✖ dplyr::filter() masks stats::filter() ✖ dplyr::lag() masks stats::lag() ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors Attaching SeuratObject Attaching package: ‘cowplot’ The following object is masked from ‘package:lubridate’: stamp Loading required package: grid ======================================== ComplexHeatmap version 2.14.0 Bioconductor page: http://bioconductor.org/packages/ComplexHeatmap/ Github page: https://github.com/jokergoo/ComplexHeatmap Documentation: http://jokergoo.github.io/ComplexHeatmap-reference If you use it in published research, please cite either one: - Gu, Z. Complex Heatmap Visualization. iMeta 2022. - Gu, Z. Complex heatmaps reveal patterns and correlations in multidimensional genomic data. Bioinformatics 2016. The new InteractiveComplexHeatmap package can directly export static complex heatmaps into an interactive Shiny app with zero effort. Have a try! This message can be suppressed by: suppressPackageStartupMessages(library(ComplexHeatmap)) ======================================== ======================================== circlize version 0.4.15 CRAN page: https://cran.r-project.org/package=circlize Github page: https://github.com/jokergoo/circlize Documentation: https://jokergoo.github.io/circlize_book/book/ If you use it in published research, please cite: Gu, Z. circlize implements and enhances circular visualization in R. Bioinformatics 2014. This message can be suppressed by: suppressPackageStartupMessages(library(circlize)) ======================================== Loading required package: SingleCellExperiment Loading required package: SummarizedExperiment Loading required package: MatrixGenerics Loading required package: matrixStats Attaching package: ‘matrixStats’ The following object is masked from ‘package:dplyr’: count Attaching package: ‘MatrixGenerics’ The following objects are masked from ‘package:matrixStats’: colAlls, colAnyNAs, colAnys, colAvgsPerRowSet, colCollapse, colCounts, colCummaxs, colCummins, colCumprods, colCumsums, colDiffs, colIQRDiffs, colIQRs, colLogSumExps, colMadDiffs, colMads, colMaxs, colMeans2, colMedians, colMins, colOrderStats, colProds, colQuantiles, colRanges, colRanks, colSdDiffs, colSds, colSums2, colTabulates, colVarDiffs, colVars, colWeightedMads, colWeightedMeans, colWeightedMedians, colWeightedSds, colWeightedVars, rowAlls, rowAnyNAs, rowAnys, rowAvgsPerColSet, rowCollapse, rowCounts, rowCummaxs, rowCummins, rowCumprods, rowCumsums, rowDiffs, rowIQRDiffs, rowIQRs, rowLogSumExps, rowMadDiffs, rowMads, rowMaxs, rowMeans2, rowMedians, rowMins, rowOrderStats, rowProds, rowQuantiles, rowRanges, rowRanks, rowSdDiffs, rowSds, rowSums2, rowTabulates, rowVarDiffs, rowVars, rowWeightedMads, rowWeightedMeans, rowWeightedMedians, rowWeightedSds, rowWeightedVars Loading required package: GenomicRanges Loading required package: stats4 Loading required package: BiocGenerics Attaching package: ‘BiocGenerics’ The following object is masked from ‘package:limma’: plotMA The following objects are masked from ‘package:lubridate’: intersect, setdiff, union The following objects are masked from ‘package:dplyr’: combine, intersect, setdiff, union The following objects are masked from ‘package:stats’: IQR, mad, sd, var, xtabs The following objects are masked from ‘package:base’: anyDuplicated, aperm, append, as.data.frame, basename, cbind, colnames, dirname, do.call, duplicated, eval, evalq, Filter, Find, get, grep, grepl, intersect, is.unsorted, lapply, Map, mapply, match, mget, order, paste, pmax, pmax.int, pmin, pmin.int, Position, rank, rbind, Reduce, rownames, sapply, setdiff, sort, table, tapply, union, unique, unsplit, which.max, which.min Loading required package: S4Vectors Attaching package: ‘S4Vectors’ The following objects are masked from ‘package:lubridate’: second, second<- The following objects are masked from ‘package:dplyr’: first, rename The following object is masked from ‘package:tidyr’: expand The following objects are masked from ‘package:base’: expand.grid, I, unname Loading required package: IRanges Attaching package: ‘IRanges’ The following object is masked from ‘package:lubridate’: %within% The following objects are masked from ‘package:dplyr’: collapse, desc, slice The following object is masked from ‘package:purrr’: reduce Loading required package: GenomeInfoDb Loading required package: Biobase Welcome to Bioconductor Vignettes contain introductory material; view with 'browseVignettes()'. To cite Bioconductor, see 'citation("Biobase")', and for packages 'citation("pkgname")'. Attaching package: ‘Biobase’ The following object is masked from ‘package:MatrixGenerics’: rowMedians The following objects are masked from ‘package:matrixStats’: anyMissing, rowMedians Attaching package: ‘SummarizedExperiment’ The following object is masked from ‘package:SeuratObject’: Assays The following object is masked from ‘package:Seurat’: Assays Loading required package: scuttle
library(future)
#for 200gb ram
options(future.globals.maxSize = 200000 * 1024^2)
sessionInfo()
R version 4.2.2 (2022-10-31) Platform: x86_64-conda-linux-gnu (64-bit) Running under: CentOS Stream 8 Matrix products: default BLAS/LAPACK: /hpc/group/pbenfeylab/tmn23/miniconda3/envs/muscat/lib/libopenblasp-r0.3.21.so locale: [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8 [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8 [7] LC_PAPER=en_US.UTF-8 LC_NAME=C [9] LC_ADDRESS=C LC_TELEPHONE=C [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C attached base packages: [1] stats4 grid stats graphics grDevices utils datasets [8] methods base other attached packages: [1] future_1.31.0 scran_1.26.0 [3] scuttle_1.8.0 SingleCellExperiment_1.20.0 [5] SummarizedExperiment_1.28.0 Biobase_2.58.0 [7] GenomicRanges_1.50.0 GenomeInfoDb_1.34.8 [9] IRanges_2.32.0 S4Vectors_0.36.0 [11] BiocGenerics_0.44.0 MatrixGenerics_1.10.0 [13] matrixStats_0.63.0 limma_3.54.0 [15] muscat_1.12.0 ggrepel_0.9.3 [17] gprofiler2_0.2.1 GeneOverlap_1.34.0 [19] circlize_0.4.15 ComplexHeatmap_2.14.0 [21] cowplot_1.1.1 SeuratObject_4.1.3 [23] Seurat_4.3.0 lubridate_1.9.2 [25] forcats_1.0.0 stringr_1.5.0 [27] dplyr_1.1.0 purrr_1.0.1 [29] readr_2.1.4 tidyr_1.3.0 [31] tibble_3.1.8 ggplot2_3.4.1 [33] tidyverse_2.0.0 loaded via a namespace (and not attached): [1] pbdZMQ_0.3-9 scattermore_0.8 [3] bit64_4.0.5 irlba_2.3.5.1 [5] DelayedArray_0.24.0 data.table_1.14.8 [7] KEGGREST_1.38.0 RCurl_1.98-1.10 [9] doParallel_1.0.17 generics_0.1.3 [11] ScaledMatrix_1.6.0 RhpcBLASctl_0.23-42 [13] RSQLite_2.2.20 RANN_2.6.1 [15] bit_4.0.5 tzdb_0.3.0 [17] spatstat.data_3.0-0 httpuv_1.6.9 [19] viridis_0.6.2 hms_1.1.2 [21] evaluate_0.20 promises_1.2.0.1 [23] fansi_1.0.4 progress_1.2.2 [25] caTools_1.18.2 igraph_1.3.5 [27] DBI_1.1.3 geneplotter_1.76.0 [29] htmlwidgets_1.6.1 spatstat.geom_3.0-6 [31] ellipsis_0.3.2 backports_1.4.1 [33] annotate_1.76.0 aod_1.3.2 [35] deldir_1.0-6 sparseMatrixStats_1.10.0 [37] vctrs_0.5.2 ROCR_1.0-11 [39] abind_1.4-5 cachem_1.0.6 [41] withr_2.5.0 progressr_0.13.0 [43] sctransform_0.3.5 prettyunits_1.1.1 [45] goftest_1.2-3 cluster_2.1.4 [47] IRdisplay_1.1 lazyeval_0.2.2 [49] crayon_1.5.2 genefilter_1.80.0 [51] spatstat.explore_3.0-6 edgeR_3.40.0 [53] pkgconfig_2.0.3 nlme_3.1-162 [55] vipor_0.4.5 blme_1.0-5 [57] rlang_1.0.6 globals_0.16.2 [59] lifecycle_1.0.3 miniUI_0.1.1.1 [61] rsvd_1.0.5 polyclip_1.10-4 [63] lmtest_0.9-40 Matrix_1.5-3 [65] IRkernel_1.3.2 boot_1.3-28.1 [67] zoo_1.8-11 base64enc_0.1-3 [69] beeswarm_0.4.0 ggridges_0.5.4 [71] GlobalOptions_0.1.2 png_0.1-8 [73] viridisLite_0.4.1 rjson_0.2.21 [75] bitops_1.0-7 KernSmooth_2.23-20 [77] Biostrings_2.66.0 blob_1.2.3 [79] DelayedMatrixStats_1.20.0 shape_1.4.6 [81] parallelly_1.34.0 spatstat.random_3.1-3 [83] beachmat_2.14.0 scales_1.2.1 [85] memoise_2.0.1 magrittr_2.0.3 [87] plyr_1.8.8 ica_1.0-3 [89] gplots_3.1.3 zlibbioc_1.44.0 [91] compiler_4.2.2 dqrng_0.3.0 [93] RColorBrewer_1.1-3 clue_0.3-64 [95] lme4_1.1-31 DESeq2_1.38.0 [97] fitdistrplus_1.1-8 cli_3.6.0 [99] XVector_0.38.0 lmerTest_3.1-3 [101] listenv_0.9.0 patchwork_1.1.2 [103] pbapply_1.7-0 TMB_1.9.2 [105] MASS_7.3-58.2 tidyselect_1.2.0 [107] stringi_1.7.12 BiocSingular_1.14.0 [109] locfit_1.5-9.7 tools_4.2.2 [111] timechange_0.2.0 future.apply_1.10.0 [113] parallel_4.2.2 uuid_1.1-0 [115] bluster_1.8.0 foreach_1.5.2 [117] metapod_1.6.0 gridExtra_2.3 [119] Rtsne_0.16 digest_0.6.31 [121] shiny_1.7.4 Rcpp_1.0.10 [123] broom_1.0.3 later_1.3.0 [125] RcppAnnoy_0.0.20 httr_1.4.4 [127] AnnotationDbi_1.60.0 Rdpack_2.4 [129] colorspace_2.1-0 XML_3.99-0.13 [131] tensor_1.5 reticulate_1.28 [133] splines_4.2.2 statmod_1.5.0 [135] uwot_0.1.14 spatstat.utils_3.0-1 [137] scater_1.26.0 sp_1.6-0 [139] plotly_4.10.1 xtable_1.8-4 [141] jsonlite_1.8.4 nloptr_2.0.3 [143] R6_2.5.1 pillar_1.8.1 [145] htmltools_0.5.4 mime_0.12 [147] glue_1.6.2 fastmap_1.1.0 [149] minqa_1.2.5 BiocParallel_1.32.5 [151] BiocNeighbors_1.16.0 codetools_0.2-19 [153] utf8_1.2.3 lattice_0.20-45 [155] spatstat.sparse_3.0-0 numDeriv_2016.8-1.1 [157] pbkrtest_0.5.2 ggbeeswarm_0.7.1 [159] leiden_0.4.3 gtools_3.9.4 [161] survival_3.5-3 glmmTMB_1.1.5 [163] repr_1.1.6 munsell_0.5.0 [165] GetoptLong_1.0.5 GenomeInfoDbData_1.2.9 [167] iterators_1.0.14 variancePartition_1.28.0 [169] reshape2_1.4.4 gtable_0.3.1 [171] rbibutils_2.2.13
rc.integrated <- readRDS("../../CheWei/scRNA-seq/Integrated_Objects/rc.integrated_11S_CVP_BRI1_seu3_20230315.rds")
rc.integrated
An object of class Seurat 71788 features across 72381 samples within 3 assays Active assay: integrated (17520 features, 17520 variable features) 2 other assays present: RNA, SCT 4 dimensional reductions calculated: pca, umap, umap_3D, umap_2D
# remove sc_12 dc1 dc2
rc.integrated <- subset(rc.integrated, subset = orig.ident %in% c("briT", "briTR","sc_130", "sc_131", "sc_132", "sc_134", "sc_135", "sc_136"))
table(rc.integrated$orig.ident)
briT briTR sc_130 sc_131 sc_132 sc_134 sc_135 sc_136 7483 8318 6589 7615 6550 7745 5038 6089
rc.integrated$geno <- rc.integrated$orig.ident
rc.integrated$geno <- gsub("briT","bri1_T",rc.integrated$geno)
rc.integrated$geno <- gsub("bri1_TR","pCVP2_BRI1_Citrine_bri1_T",rc.integrated$geno)
rc.integrated$geno <- gsub("sc_130","WT",rc.integrated$geno)
rc.integrated$geno <- gsub("sc_131","bri1_T",rc.integrated$geno)
rc.integrated$geno <- gsub("sc_132","pCVP2_BRI1_Citrine_bri1_T",rc.integrated$geno)
rc.integrated$geno <- gsub("sc_134","WT",rc.integrated$geno)
rc.integrated$geno <- gsub("sc_135","bri1_T",rc.integrated$geno)
rc.integrated$geno <- gsub("sc_136","pCVP2_BRI1_Citrine_bri1_T",rc.integrated$geno)
rc.integrated$source <- rc.integrated$orig.ident
rc.integrated$source <- gsub("briT","Graeff-2021",rc.integrated$source)
rc.integrated$source <- gsub("bri1_TR","Graeff-2021",rc.integrated$source)
rc.integrated$source <- gsub("Graeff-2021R","Graeff-2021",rc.integrated$source)
rc.integrated$source <- gsub("sc_130","Benfey",rc.integrated$source)
rc.integrated$source <- gsub("sc_131","Benfey",rc.integrated$source)
rc.integrated$source <- gsub("sc_132","Benfey",rc.integrated$source)
rc.integrated$source <- gsub("sc_134","Benfey",rc.integrated$source)
rc.integrated$source <- gsub("sc_135","Benfey",rc.integrated$source)
rc.integrated$source <- gsub("sc_136","Benfey",rc.integrated$source)
rc.integrated$rep <- rc.integrated$orig.ident
rc.integrated$rep <- gsub("briT","1",rc.integrated$rep)
rc.integrated$rep <- gsub("bri1_TR","1",rc.integrated$rep)
rc.integrated$rep <- gsub("sc_130","2",rc.integrated$rep)
rc.integrated$rep <- gsub("sc_131","2",rc.integrated$rep)
rc.integrated$rep <- gsub("sc_132","2",rc.integrated$rep)
rc.integrated$rep <- gsub("sc_134","3",rc.integrated$rep)
rc.integrated$rep <- gsub("sc_135","3",rc.integrated$rep)
rc.integrated$rep <- gsub("sc_136","3",rc.integrated$rep)
table(rc.integrated$orig.ident, rc.integrated$geno)
bri1_T pCVP2_BRI1_Citrine_bri1_T WT
briT 7483 0 0
briTR 0 8318 0
sc_130 0 0 6589
sc_131 7615 0 0
sc_132 0 6550 0
sc_134 0 0 7745
sc_135 5038 0 0
sc_136 0 6089 0
table(rc.integrated$orig.ident, rc.integrated$source)
Benfey Graeff-2021
briT 0 7483
briTR 0 8318
sc_130 6589 0
sc_131 7615 0
sc_132 6550 0
sc_134 7745 0
sc_135 5038 0
sc_136 6089 0
feature_names <- read_tsv("../data/features.tsv.gz", col_names = c("AGI", "Name", "Type")) %>%
select(-Type) %>%
distinct()
Rows: 32833 Columns: 3 ── Column specification ──────────────────────────────────────────────────────── Delimiter: "\t" chr (3): AGI, Name, Type ℹ Use `spec()` to retrieve the full column specification for this data. ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
rc.integrated$geno <- factor(rc.integrated$geno, levels=c("WT", "bri1_T", "pCVP2_BRI1_Citrine_bri1_T"))
time_zonecell_typetime_zone_cell_typetime_zone_cell_subtypesrc.integrated$cell_type <- rc.integrated$celltype.anno.Li.crude
rc.integrated$time_zone <- rc.integrated$time.anno
rc.integrated$time_zone_cell_type <- rc.integrated$time.celltype.anno.Li.crude
table(rc.integrated$orig.ident, rc.integrated$cell_type)
Atrichoblast Columella Cortex Endodermis Lateral Root Cap Pericycle
briT 640 403 332 401 2368 1750
briTR 700 637 228 442 4148 1167
sc_130 1445 895 972 725 1058 143
sc_131 1441 672 777 490 1800 812
sc_132 1089 474 720 544 1597 554
sc_134 1171 582 906 573 2101 633
sc_135 1055 294 544 480 991 428
sc_136 1176 649 878 617 1240 247
Phloem Procambium Quiescent Center Trichoblast Xylem
briT 204 815 7 200 363
briTR 227 358 13 148 250
sc_130 52 127 0 1000 172
sc_131 130 320 0 900 273
sc_132 138 288 0 916 230
sc_134 123 348 0 1041 267
sc_135 88 260 0 707 191
sc_136 78 186 0 871 147
table(rc.integrated$time.anno)
Distal Columella Distal Lateral Root Cap Elongation
3943 5190 15349
Maturation Meristem Proximal Columella
10718 9559 633
Proximal Lateral Root Cap
10035
# Plot celltype annotation Li
order <- c("Quiescent Center", "Ground Tissue","Columella", "Lateral Root Cap", "Atrichoblast", "Trichoblast", "Cortex", "Endodermis", "Phloem","Protophloem", "Xylem", "Procambium","Pericycle","Phloem Pole Pericycle", "Protoxylem", "Metaxylem", "Unknown")
palette <- c("#9400D3", "#DCD0FF","#5AB953", "#BFEF45", "#008080", "#21B6A8", "#82B6FF", "#0000FF","#E6194B", "#DD77EC", "#9A6324", "#FFE119", "#FF9900", "#FFD4E3", "#9A6324", "#DDAA6F", "#EEEEEE")
rc.integrated$cell_type <- factor(rc.integrated$cell_type, levels = order[sort(match(unique(rc.integrated$cell_type),order))])
color <- palette[sort(match(unique(rc.integrated$cell_type),order))]
options(repr.plot.width=15, repr.plot.height=7)
(BR_cell <- DimPlot(rc.integrated, reduction = "umap", group.by = "cell_type", cols = color, split.by = 'geno', pt.size = 0.75, ncol=3))
options(repr.plot.width=15, repr.plot.height=7)
FeaturePlot(rc.integrated, features="BRI1-mCitrine", split.by = "geno", order=T)
Warning message in FeaturePlot(rc.integrated, features = "BRI1-mCitrine", split.by = "geno", : “All cells have the same value (0) of BRI1-mCitrine.”
options(repr.plot.width=30, repr.plot.height=7)
(BR_cell <- DimPlot(rc.integrated, reduction = "umap", group.by = "cell_type", cols = color, split.by = 'orig.ident', pt.size = 0.75, ncol=8))
DefaultAssay(rc.integrated) <- "SCT"
FeaturePlot(rc.integrated, features="BRI1-mCitrine", split.by = "orig.ident", order=T)
Warning message in FeaturePlot(rc.integrated, features = "BRI1-mCitrine", split.by = "orig.ident", : “All cells have the same value (0) of BRI1-mCitrine.” Warning message in FeaturePlot(rc.integrated, features = "BRI1-mCitrine", split.by = "orig.ident", : “All cells have the same value (0) of BRI1-mCitrine.”
saveRDS(rc.integrated, file = "/hpc/group/pbenfeylab/CheWei/scRNA-seq/Integrated_Objects/rc.integrated_8S_CVP_BRI1_seu3_annotated_20230316.rds")